############# Market Share Models High level FOR within funders

rm(list=ls())

library(here)
here::here()  # points to project root

library(tidyverse)
`%nin%` <- Negate(`%in%`)

data <- read_csv("us_funding_by_funder_by_year.csv")

for_grant <- read_csv("us_funding_by_funder_by_year_high_level_for_field.csv")
for_grant <- for_grant %>% filter(funder %nin% c("California Institute for Regenerative Medicine"))



for_grant_funding <- for_grant %>% group_by(for_field, funder, start_year) %>% dplyr::summarise(funding_usd_adj = sum(funding_usd_adj))

for_grant_covars <- for_grant %>% dplyr::select(-`...1`,-funding_usd, -funding_usd_adj, -pct_change, -abs_change) %>% unique()

for_grant <- for_grant_funding %>% left_join(for_grant_covars)


for_grant <- for_grant %>% mutate(cong_D = case_when(majority_house == 1 & majority_sen == 1 ~ 1,
                                                     TRUE ~ 0), 
                                  trifecta_D = case_when(majority_house == 1 & majority_sen == 1 & president_party == 1 ~ 1,
                                                         TRUE ~ 0),
                                  cong_R = case_when(majority_house == 0 & majority_sen ==  0 ~ 1,
                                                     TRUE ~ 0), 
                                  trifecta_R = case_when(majority_house == 0 & majority_sen == 0 & president_party == 0 ~ 1,
                                                         TRUE ~ 0))

library("Hotelling")
grant_years <- for_grant %>% group_by(start_year, funder) %>% dplyr::summarise(yearly_funding = sum(funding_usd_adj))
for_grant <- for_grant %>% left_join(grant_years)
for_grant$funding_share <- (for_grant$funding_usd_adj+1)/(for_grant$yearly_funding+1)


for_grant_rank_dir_DV <- for_grant %>% dplyr::select(for_field, funder, start_year,funding_share,trifecta_D, trifecta_R,cong_R, cong_D)


for_grant_rank_dir_DV <- for_grant_rank_dir_DV %>% filter(start_year > 1994) %>% ungroup() %>% unique() %>% pivot_wider(id_cols = c(start_year,funder, trifecta_D, trifecta_R,cong_R, cong_D), names_from = for_field,  values_from = funding_share)

for_grant_rank_dir_covars <- for_grant_rank_dir_DV[,1:6]
for_grant_rank_dir_DV <- for_grant_rank_dir_DV[,7:28]
for_grant_rank_dir_DV[is.na(for_grant_rank_dir_DV)] <- min(for_grant_rank_dir_DV[for_grant_rank_dir_DV>0])

for_grant_rank_dir_DV <- clr(for_grant_rank_dir_DV)

for_grant_rank_dir <- bind_cols(for_grant_rank_dir_DV, for_grant_rank_dir_covars)

for_grant_rank_dir <- for_grant_rank_dir %>% pivot_longer(cols = 1:22, names_to = "FOR_text", values_to = "CLR_Share")

for_grant_rank_dir <- for_grant_rank_dir %>% mutate(cong_control = case_when(cong_D == 1 ~ "Dem",
                                                                             cong_R == 1 ~ "Rep",
                                                                             T ~ "Divided"))

for_grant_rank_dir <- for_grant_rank_dir %>% mutate(total_control = case_when(trifecta_D == 1 ~ "Dem",
                                                                              trifecta_R == 1 ~ "Rep",
                                                                              T ~ "Divided"))

for_grant_rank_dir$cong_control <- as.factor(for_grant_rank_dir$cong_control)
for_grant_rank_dir$total_control <- as.factor(for_grant_rank_dir$total_control)


market_share_m1 <- lm(CLR_Share ~  cong_control*FOR_text, data=for_grant_rank_dir)
library(interactions)
msp1 <- cat_plot(market_share_m1, pred = "FOR_text", modx ="cong_control") + coord_flip() + ggtitle("Field of Research Market Share within Funder") +
  theme(legend.position = "Bottom")

msp1$data %>% ggplot(aes(y=fct_reorder(FOR_text,CLR_Share, "median" ), x= CLR_Share, xmin = ymin, xmax = ymax, group = cong_control, color = cong_control)) + geom_pointrange(position = position_dodge(width = .7), alpha = .5) + 
  scale_color_manual("Party Control of Congress",  values = c("#156B90","#666666", "#9A3E25"))   + theme_minimal() + theme(legend.position = "bottom") +
  xlab("Centered Log Ratio\n(Within Funder)") + ylab("") -> p1


############# Market Share Models High level FOR aggregated across funders
for_grant <- read_csv("us_funding_by_funder_by_year_high_level_for_field.csv")
for_grant <- for_grant %>% filter(funder %nin% c("California Institute for Regenerative Medicine"))



for_grant_funding <- for_grant %>% group_by(for_field, start_year) %>% dplyr::summarise(funding_usd_adj = sum(funding_usd_adj))

for_grant_covars <- for_grant %>% dplyr::select(-`...1`,-funder,-funding_usd, -funding_usd_adj, -pct_change, -abs_change) %>% unique()

for_grant <- for_grant_funding %>% left_join(for_grant_covars)

for_grant_rank <- for_grant %>% group_by(start_year) %>% mutate(rank = min_rank(desc(funding_usd_adj))) 

grant_years <- for_grant_rank %>% group_by(start_year) %>% dplyr::summarise(yearly_funding = sum(funding_usd_adj))
for_grant_rank <- for_grant_rank %>% left_join(grant_years)
for_grant_rank$funding_share <- for_grant_rank$funding_usd_adj/for_grant_rank$yearly_funding

for_grant_rank <- for_grant_rank %>% mutate(cong_D = case_when(majority_house == 1 & majority_sen == 1 ~ 1,
                                                               TRUE ~ 0), 
                                            trifecta_D = case_when(majority_house == 1 & majority_sen == 1 & president_party == 1 ~ 1,
                                                                   TRUE ~ 0),
                                            cong_R = case_when(majority_house == 0 & majority_sen ==  0 ~ 1,
                                                               TRUE ~ 0), 
                                            trifecta_R = case_when(majority_house == 0 & majority_sen == 0 & president_party == 0 ~ 1,
                                                                   TRUE ~ 0))
for_grant_rank <- for_grant_rank %>% mutate(cong_control = case_when(cong_D == 1 ~ "Dem",
                                           cong_R == 1 ~ "Rep",
                                           T ~ "Divided"))

for_grant_rank %>% filter(start_year > 1994) %>% dplyr::select(for_field, cong_control, start_year,funding_usd_adj ) -> field_party_year

for_grant_rank %>% filter(start_year > 1994) %>% dplyr::group_by(start_year) %>% dplyr::summarise(tot_fund = sum(funding_usd_adj, na.rm = T)) -> year_sum

field_party_share <- field_party_year %>% left_join(year_sum) %>% mutate(fund_share = funding_usd_adj/tot_fund) %>% group_by(for_field, cong_control) %>%dplyr::summarise(mean_share = mean(fund_share))

ggplot(field_party_share, aes(x=mean_share, y=fct_reorder(for_field, mean_share, "median") , group = cong_control, fill = cong_control)) + geom_bar(stat="identity", position = "dodge") + theme_minimal() + 
  theme(legend.position = "bottom") + scale_fill_manual("Party Control of Congress",  values = c("#156B90","#666666", "#9A3E25")) +
  xlab("Average Share of Yearly Grant Outlays") + ylab("") -> p3


ggsave("Fig5.pdf", p3, width = 7, height = 8)

cowplot::plot_grid(p1,p2, nrow = 1) -> out_plot

ggsave("FigS5BC.pdf", out_plot, device = "pdf", width = 12, height = 6, units="in")



